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ABSTRACT 

Pulsar timing observations are used to place constraints on the rate of coalescence of su- 
permassive black-hole (SMBH) binaries as a function of mass and redshift. In contrast to the 
indirect constraints obtained from other techniques, pulsar timing observations provide a direct 
constraint on the black-hole merger rate. This is possible since pulsar timing is sensitive to the 
gravitational waves (GWs) emitted by these sources in the final stages of their evolution. We 
find that upper bounds calculated from the recently published Parkes Pulsar Timing Array data 
are just above theoretical predictions for redshifts below 10. In the future, with improved timing 
precision and longer data spans, we show that a non-detection of GWs will rule out some of the 
available parameter space in a particular class of SMBH binary merger models. We also show 
that if we can time a set of pulsars to 10 ns timing accuracy, for example, using the proposed 
Square Kilometre Array, it should be possible to detect one or more individual SMBH binary 
systems. 

Subject headings: pulsars: general — gravitational waves — methods: data analysis — early 
universe — Galaxies: statistics 



1. Introduction 



Pulsar timing observations (for a review of the techniques, see Lorimer & Krame JboOS : Edwards et al. 



2006h provide a unique opportunity to study low-frequency (10 ^-10 ^ Hz) gravitational waves (GWs; 
Sazhin"l978':'Detweileij|l979l: lBertotti et al.lll983l:lHem ngs & Downs'l983':'Foster & B acker" 1 990*: 'Kas pi et al 



1994 : .Jenet et al...2005,) . Previous work dRomani & Tavlor l983 : .Kaspi et al.J994. : .Lommen2002. : Jenet et al 
2006]) placed upper limits on a stochastic background of GWs. These limits were reported in terms of ei- 



ther the amplitude of the GW characteristic strain spectrum, hc{f), or the normalized GW energy density, 

^gw(/)- 
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In recent years, researchers have proposed that supermassive black-hole (SMBH) binary systems dis- 



Jaffe & Backer 


2003; 


Wvithe & Loeb 


2003; 


Enoki et al. 


2004; 


Sesana et al. 


2004 



or 



non-detection, of such GWs provides a constraint on the rate of coalescence of SMBH binary systems. We 
emphasize that such constraints are model-independent as opposed to the indirect constraints that can be 
inferred from observed galaxy distributions. 

We will show that existing pulsar data sets do not provide stringent constraints on the coalescence 



rate. H owever, future data sets from the Parkes Puls ar Timing Array (P PTA; lManchestetil2008l;lHobbs et al. 



2009al) and similar North American (NANOGrav; Lfenet et all 120091) and European timing array (EPTA; 
Stappers et al.ll2006r) projects aim to produce data sets on 20 or more pulsars with root-mean-square (rms) 
timing residuals close to 100 ns. In the longer term, we expect that pulsar timing array projects using future 
telescopes such as the Square Kilometre Array (SkaA will be able to time many hundreds of pulsars with 
exquisite timing precision. 

The outline of this paper is as follows. In § 2, the physics of GW emission from SMBH binary systems 
is reviewed together with the effects of GWs on pulsar timing. We describe how to constrain the coalescence 
rate using two different techniques, one valid when there is a large number of expected sources and the other 
valid when only a few sources are expected. In § 3 we show the recent and projected rate constraints for 
various different observing systems. These observationally constrained rates are then compared to the rates 
implied by local galaxy-merger observations. This work is summarized in section § 4. 



2. Constraining the SMBH Merger Rate 



We define a SMBH as a black h ole with mass greater than 10^ M(T). There is abund ant evidence that 
such SMBHs exis t both nearby (e .g., iKormendy & Richstondll992l ; iMiyoshi et al.lll995f) and at high red- 
shifts z ~ 6 (e.g.. Fan et al.ll200lh . Two orbiting SMBHs resulting from a galaxy merger would emit large 
amplitude GWs. Such GW sources are important targets for spa ce-based detectors such as the Laser In- 
terferometer Space Antenna (LISA) and pulsar timing arrays (e.g. JWyithe & Loebll2003l ; Enoki et al.ll2004l ; 
Sesana et al.ll2004r) . Observational GW astronomy can be used to resolve whether SMBH binary systems 
form and whether the two black holes can get close enough to emit detectable GWs. 

This work focuses on the rate constraints measurable by pulsar timing observations. In order to deter- 
mine the co alescence rate , one needs to know the expected GW amplitude emitted by a SMBH binary. This 
is given by (|Thomdll987h 



[7Tf{l + z)Y 



(1) 



5 c^D{z) 

where is the "chirp mass" of the SMBH binary given by = (MiM2)^/^(Afi + Ms)"^/^, Mi and 
M2 are the individual black hole masses, / is observed GW frequency, D{z) is the comoving distance to 



'See www.skatelescope.org 
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the system 

c 



with E{z) = y^r^A + ^m(l + for a ACDM cosmofogical model (hereafter we adopt Hq = 72 km s ^ Mpc 
Jim = 0.3, JIa = 0.7). GWs from such a source w ill induce sinusoidally oscillating arrival-time variations 



whose amplitude, At, is given by (IJenet et al.ll2004 and references therein): 



At = — [1 + cos(6')] sin(2(/.) siu{ujD[l - cos(6')]/2c}, (3) 

U) 

where oj is the GW frequency in radians/s, 9 is the angle on the sky between the pulsar direction and the GW 
source direction, is the GW polarization angle, and D is the distance to the pulsar. The maximum induced 
timing residual amplitude, /is/o;, is plotted in Figure [Has a function of redshift for chirp masses of 10^ and 
10^" Mq and observed GW frequencies of (1 year)~^ and (10 year)"^. The most notable feature in this 
Figure is that these amplitude curves are not monotonically decreasing with increasing redshift. The reason 
for this is that the observing frequency is held fixed and the frequency in the frame of the emitting system 
increases with increasing redshift. For a binary system, the emitted GW amplitude increases with increasing 
frequency. For large z, this increase is faster than the decrease due to increasing comoving distance D{z). 
Also note that the curves cutoff at large z. This is because there is a maximum orbital frequency allowed 
before the blac k holes plunge together. This maximum frequency was taken to be c^/ (6^/^7rGM) assuming 



a circular orbit (lHughesll2002l) . Here, M = Mi + M2 is the total mass of a binary system. 



Rms timing residuals (for a typical 1-hour observation) for the best pulsar data sets are currently around 
100 ns. Multiple observations combined with improved systems will bring the effective sensitivity down to 
around 10 ns. In this case. Figure [U shows that pulsar timing will be sensitive to individual SMBH binary 
systems with chirp masses greater than about 10^ Mq. Figure [U also shows that this sensitivity extends 
to large redshifts. This fact greatly increases the chances of detecting individual sources. An ensemble of 
lower-mass SMBH binary systems will be detectable as a stochastic background if there is a large enough 
population of sources. 

Since pulsar timing techniques are sensitive to SMBH binary systems up to high redshifts, the non- 
detection of any sources, either individual binaries or a background generated by an ensemble of binary 
systems, may be used to place a direct constraint on (PR/dMcdz, the sky-averaged rate of coalescence of 
binary SMBHs per unit chirp mass, Mc, per unit redshift z. The total number of binary SMBHs with chirp 
mass between Mc and + AMc and located between z and z + Az merging between time t and t + At 
is given by AMcAzAtcfR/dMcdz. Constraints placed on this quantity may be used to rule out various 
binary SMBH formation models. 

We present two methods for determining the coalescence rate from pulsar timing data. The first is 
valid when the rate is high enough that the GWs form a stochastic background (the "stochastic constraint") 
and there is a large number of sources per resolvable frequency bin. The second method (the "Poissonian 
constraint") provides an estimate of the coalescence rate when the stochastic constraint does not hold. For a 
real data set it is practical first to assume the stochastic constraint, determine the coalescence rate and check 
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whether the rate is high enough for the assumption to be valid. If not then the Poissonian constraint should 
be used. 



2.1. Stochastic constraint 



Here, it is assumed that a large number of SMBH bin ary sources form an in coherent background of 
GWs. The power spectrum of such a background is given by Uaffe & Backeii (l2003h 



pii) 



-1 



dzclMc, 



(4) 



where hs{f, M^, z) is given by Equation [T] and df /dt is the rate of change of the obser ved GW frequency. 
For th e case of a binary system evolving under general relativity alone, this is given by (IPeters & Mathews 
1963h 

df 96(GM^y (^^)8/3 ^(1^^)5/3 



(5) 



dt 5 V c3 

Typically, a stochastic background of GWs is described by its characteristic strain spectrum, hc{f), which 
is assumed to take on a power law form: 



(6) 



where /yr = 1/ (lyear) and A is the characteristic strain at a period of one y ear. For a background generated 
by SMBH binaries, a = -2/3 in frequenc y range lO^^-lO^'^ Hz (..Taffe & Backer! boosi : Iwvithe & Loeb 



20031 : ISesana et al.ll2004l : 



Enoki et al. 



20041) . Note that the characteristic strain spectrum is related to the 



power spectrum by P{f) = hcifY/f- 

Pulsar tuning data sets provide an upper bound, ^up> on A. Such bounds limit the power spectrum of 
the GW strain. The upper bound, -Pup(/), may be written as 



42 



/yr ^/yr' 



(V) 



Since this is an upper bound, we have 

P{f) < Pupif). (8) 

In order to use Equations |4] and [8] to obtain a constraint on the differential rate of coalescence itself, the 
integrand is rewritten in an equivalent form: 



P{f) 



00 roo 



hsif,M„z) 



d^R 



dlg{l + z)dlg{Mc) \dt 



df 



dlg(l + z)dlg(Me 



(9) 



Note that both d'^R/dlg{l + z)dlg(Mc) and df /dt depend on z and Mc, although the explicit dependence 
is not written. A constraint on Pif ) is a direct constraint on the integral in the above expression. In 
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order to obtain an estimate of the upper bound on the integrand, we follow the same Une of reasoning used 
to place constrai nts on the differ ential energy density of gravitational waves using bounds from big bang 
nucleosynthesis jlVIagg iore 2000 ). First, we note that the limits in the integral of equation |9]are from to cx) 
for lg(l + z) and from —oo to oo for Ig(Mc). Consider this integral over a small region bounded by lg{Mcj ) 
to lg(Mc2) and lg(l + zi) to lg(l + 22)- Denote this integral as Ps(/) where 

lg(M.,) ^lg(l+.2) ^2^ fdfY^ 



Ps{f)= / hs{f,M„zy [-^ dlg{l + z)dlg{M,). (10) 

Jig{M,^) Jig{i+z^) dlg{l + z)dlg{Mc) \dt J 

The mean value theorem tells us that there exist values of and z, written as M* and z*, such that 

Ps{f) = hs{f,M:,z*) \,^ ^'^'^^^ 'Alg(l + z)Alg(Me), (11) 



dlg{l + z)dlg{M,) ydt 

where Alg(l + z) = lg(l + Z2) - lg(l + zi) and Alg(Mc) = lg(Mc2) - lg(McJ. Next, we assume that 
the integrand varies slowly over the region of integration so that Ps{f) does not change much as long as M* 
and z* are chosen within this region. From this, we see that Ps{f) is approximately given by 

P,{f)^hs{f,M,,zf——^-^——(^) \lg{l + z)Alg{M,) (12) 



dlg{l + z)dlg{Mc) \dt 

for any value of Mc G [M^jMcj] and z G [^1,2:2]- Next, since the integrand in equation |9] is positive 
definite, it follows that 

Psif) < P{f). (13) 

From equations 171 [8l [T2l and [13] we have: 



hs{f\M,,zY (M\ \7;) Alg(l + z)Alg(Afe)<^h^ , (14) 

d\g{l + z)d\g[Mc) \dt J Jyr ^/yr^ 

where we are free to choose any value for Mc and z provided that the integrand is slowly varying over the 
appropriate region. Using the known expressions for hs{f,Mc,z) and df /dt, one can use equation [141 to 
obtain a constraint on the differential rate of coalescence. Assuming A Ig(Mc) = 1 and A lg(l + z) = 0.2, 
the constraint takes the form 



dlg(l + z)(ilg(Me) "P (GM,)5/3 

which is independent of frequency. Therefore, a measured bound, A^p, can be used directly to constrain the 
SMBH coalescence rate. Note that, if one believes that the integrand changes more rapidly with z and/or 
Mc, one can choose a sufficiently small integration range over which this assumption is true and then rescale 
the above constraint. 

This constraint is only valid when there are a large number of sources emitting into the same fre- 
quency band at the same time. In order for this to be true, the GW amplitude of each source must be much 
less than the minimum detectable amplitude as determined by the statistical properties of the pulsar tim- 
ing data, otherwise a detection would have been made. This reasoning leads to the following constraint: 
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hs{f, Mc, zf' <^ Pup(/)A/, where A/ is the resolution bandwidth which is taken to be 1/Tobs- For the 
purposes of making numerical estimates, the following constraint is used 

hsif,M„zf<0.1^^. (16) 

Jobs 

For a fixed chirp mass and frequency, this expression is a constraint on z. The most stringent constraint 
occurs when / = l/Tobs, the lowest observable frequency. Combining the above with Equations [T] and |7] 
yields the following constraint on the redshift 

The factor (1 + z)'^^^ /D{z) decreases with z until z =2.65, where it starts to increase. Hence, there is a 
bounded redshift interval over which the stochastic constraint is valid. For systems outside of this range, the 
stochastic rate limit is not valid and the Poisson rate limit discussed in the next section must be employed. 



2.2. Poisson constraint 

For this case, the sources are not numerous enough to form a stochastic background, so they must be 
treated as individual events. Assuming Poissonian statistics for the probability of an event occurring, the 
probability that no events are detected is given by e~^^^ , where (N) is the expected number of events. Since 
no events are detected in the pulsar timing data, the upper limit on the expected number, (N^^), is set so that 
— 0.05. Hence, (A^) < (A^*) = 3. If the actual expected number were greater than 3, then at least 
one source would have been detected with 95% probability. 

Provided that the expected number of events that occur within the resolution bandwidth is less than 
one, the expected number of detectable events is given by 

where Pd{Mc, z, f) is the probability of detecting a SMBH binary with chirp mass Mc at a redshift of 
z with observable frequency /. Pd takes into account non-GW noise sources that can reduce the GW 
detection efficiency. Following the same argument used in the previous section to obtain an upper bound on 
the differential rate using an integral constraint, we find that: 

< — ^.r^^ . (.9) 



As with the stochastic constraint, it is assumed that Alg(Mc) = 1 and Alg(l + z) = 0.2. This constraint 
should be appropriately rescaled if the integrand in equation [18] varies over shorter intervals. 

The above constraint requires a knowledge of Pd, the probability of det ecting a SMBH binar y system 



using pulsar timing data. The P^ is calculated using a method similar to that in lYardley et al.l (120101) . For our 
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analysis we use a Monte-Carlo simulation together with the data analysis pipeline (iHobbs et al. I l2009bl) used 



to search for GWs in real and simulated pulsar timing data. We use a Neyman-Pearson decision technique 
together with a Lomb periodogram to determine the probability of detection. For a set of Np pulsars, the 
power spectra of the timing residuals are calculated and added together. This summed power spectrum is 
used as the detector. We determine noise levels and hence detection thresholds for each frequency channel so 
that the false alarm rate for a detection is 0.001 across the entire power spectrum. The detection thresholds 
are found by producing 100,000 fake data sets by shuffling the input timing residuals, carrying out standard 
pulsar timing fits and forming the summed power spectrum. 

Once the thresholds are determined, the probability of detecting a GW with a given strain amplitude 
is calculated as follows. A GW strain amplitude A and frequency / are chosen (this procedure is repeated 
on a logarithmically spaced grid where A ranges from lO^^^-lO^^'^ and / ranges from 1/(30 years )-l/(2 
weeks)). The GW polarizatio n properties are chosen to be consistent with GWs emitted from a binary system 
and simulated as described in lHobbs et al.. ( 2009b) . The direction of the GW wave- vector is chosen from a 
distribution that is uniform on the sky while its polarization is drawn from a distribution of randomly oriented 
binary systems. The induced timing residuals from this GW source are added to a shuffled version of the 
original residuals and the summed periodogram is calculated. This is repeated 1000 times and Pd{A, f) is 
given by the number of times that the GW was detected (i.e. produced power above the threshold) divided 
by the total number of trials. 



3. Results and Discussion 

The expressions in the previous sections can be used to provide constraints on the rate of coalescence 
with any measurement of (for the stochastic case) or (for the Poissonian case). Here we discuss 



the implications of the value of presented by IJenet et al.l (l2006h . use the same data set to determine 



and simulate data sets predicting possible future timing residuals. For these future data sets we simulate (1) 
a realistic goal for existing pulsar timing array experiments (20 pulsars, timed with an rms timing residual 
of 500 ns over 10 years), (2) the goal of the PPTA project (20 pulsars, timed with an rms of 100 ns over 
5 years), and a more challenging goal of 20 pulsars, timed with an rms of 100 ns over 10 years. The rms 
timing residuals that will be achieved with future telescopes, such as the SKA, is not easy to determine. It 
may be possible to time a few pulsars with exquisite precision (with rms timing residuals of 10 ns) but other 
unmodeled noise processes may make this difficult or impossible. We therefore also simulate the following 
possible future data sets: (4) 100 pulsars timed at 100 ns over five years, (5) the same for ten years, (6) 20 
pulsars timed at 10 ns for 10 years and (7) the same for 100 pulsars. 

Figures |2] and [3] plot the stochastic constraints given by equation [15] for the different observing sce- 
narios. The horizontal error bars indicate the region of lg(l + z) over which the constraint is placed. The 
solid error bars indicate that the stochastic constraint is valid, while the dotted error bars indicate that the 
stochastic validity condition is violated over the whole range. Table[I]gives the valid redshift range for each 
data set and chirp mass. As discussed above, it was assumed that A lg(l + z) = 0.2 and A Ig(Mc) = 1 in 



0.01 0.1 1 10 100 



Redshift, z 

Fig. 1. — Induced timing residual amplitudes versus redshift of a system with a given observed frequency 
and chirp mass. 



Table 1. Upper limits on the amplitude of the stochastic GW background. A '-' indicates that there is no 

valid range for that scenario. 



Data set 




A 


up 




Valid redshift range 












(109 Mo) 


(IQio Mo) 


PPTA data Qenet et al. 2006) 


1.1 


X 


10- 


-14 


0.01-180.08 


1.78^.02 


20 PSRs-500ns-10yr 


1.1 


X 


IQ- 


-15 


0.03-602.60 




20 PSRs-lOOns-5 yr 


9.9 


X 


10- 


-16 


0.07-294.97 




20 PSRs- 100 ns- 10 yr 


2.2 


X 


10- 


-16 


0.14-115.98 




20 PSRs- 10 ns- 10 yr 


2.0 


X 


10- 


-17 






100 PSRs- 100 ns-5 yr 


5.7 


X 


10- 


-16 


0.14-121.46 




100 PSRs-lOOns-lOyr 


1.3 


X 


10- 


-16 


0.27^5.30 




100 PSRs-lOns-lOyr 


8.8 


X 


10- 


-18 
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order to calculate the constraints shown. The constraints should be rescaled if other values are assumed. 



For redshifts where the stochastic technique is not valid, the Poissonian constraint may be used. Using 
Pd{A, /), Equation[T]is used to calculate Pd{Mc, z, /) and then the right hand side of Equation [T9] is eval- 
uated numerically to determine the constraint on the differential coalescence rate. The results are plotted 
versus redshift for different chirp masses in Figures Hand |5] We note that the recently published PPTA data 
set, that with 20 pulsars timed with an rms of 500 ns over 10 years, and that with 20 pulsars timed with an 
rms of 100 ns over 5 years are not constraining for SMBH binaries with Mc = 10^ Mq and are therefore 
not plotted in the left-hand panel of Figure ID 

In order to understand how constraining the pulsar rate limits are, plots of the expected SMBH binary 
coalescence rate are also shown in the figures. Several authors have developed analytical and numerical 
techniques to estimate the coalescence rate of SMBH binaries. Here, we compare the measured rate con- 
straints to the rates predicted by the models of Uaffe & Backen (l2003b and ISesana et all (120081 120091) . For 
the case of 



Jaffe & Backen (l2003h . the differential rate of SMBH binary coalescence is given by 



47rc D{z)^ $(lg(Me)) 



dlg(l + z)dlg(Me) Ho\g{e)E{z) 



n 



(20) 



where is the total merger rate of SMBH binaries per unit comoving volume, $(lg(Mc)) is the chirp 
mass distribution of merging SMBH binaries and n is the number density of SMBH binaries given by 
n = f <I>(lg(Mc))(ilg(Mc). It is assumed that the merger rate of SMBH binaries is given by a fraction, e, 
of the galaxy merger rate and that the rate evolves as a power of (1 + z). Hence, one can write = 
eKc,(0)(l + z)^ where Kg(0) is the local merger rate of galaxy pairs and 7 is the evolution index which is 
thought to be w it hin the range — 1 < 7 < 3 (e.g 



Patton et al. 



2OO2I : iLin et al.ll2004l : iKartaltepe et al.ll2007l : 



Lin et al.ll2008l ). Iwen et al. (2009) determined 'Sig(O) and ^flgf M^)) for luminous galaxies by analyzing 
data from the Sloan Digital Sky Survey (SDSSi lYork et aLlboOoll . They found that 5Rg(0) = (1.0 ± 0.4) x 

IQ-^ Mpc-3 Gyr-i and 



lg[$(lg(M,))] = (21.7 ± 4.2) - (3.0 ± 0.5) Ig {MJMq) . 



(21) 



Note that 



Wen et all (|2009h showed that the rate implied by equation [20] together with the above estimates 



for Rg{0) and $(lg(Mr)) yield an expected characteristic strain spectrum consistent with oth er published 
estimates of he dwyithe & Loebll2003l : [sesana et al.lEo04l : lEnoki et al.ll2004l : Isesana et allboosl) . 



For the lSesana et all (|2008L l2009h models, the rate may be estimated from the following expression: 



dN„ 



1 diV (1 + z) 



dlg{l + z)dlg{Me) d\g{Me) N dz lg(e) 



(22) 



where dNm/dlg{Mc) is the mass function of coalescing SMBHs in the notation o flSesana et all (120091) . and 
dN /dz is the SMBH binary coalescence rate per unit redshift in the notation of ISesana et all (120081) . The 
constant N is given by 

'■^ dN 



N 



dz 



-dz. 



(23) 
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The limits of the integral are set by the data presented in figure 12 of ISesana et all (l2008b. The f a ctor o f 
{1 + z)/ lg(e) is used to convert the differential dz into d\g{l + z). From figure 1 of ISesana et all (120091) . 
we can estimate the maximum and minimum predicted values of the mass function dNm/dlg{Mc) over the 
four models presented therein. We find that, for Mc = 10^ Mq, the mass function lies between 10~^ yr~^ 
and 6 x 10~ ^ yr~\ - These correspond to the "Tr-SA" and the "La-SA" models, respectively, as discussed in 



Sesana et all (|2009h . For the case of Mc = W^^Mq, the predicted range lies between an 3 x 10 yr 



These values also correspond to the "Tr-SA" and "La-SA" models, respectively. F or the SMBH binary 
coalescence rate per unit redshift, dN /dz, we used the data presented in Figure 12 of ISesana et al.l (120081) . 
The three possible models shown are all approximately within a factor of two of each other. For definiteness, 
we chose the prediction based on the BVRhf model. In this case, N w 0.05 yr~^. 

Figures |2] and [3] plot the pulsar timing stochastic constraint together with the expected rates for the 
theoretical models considered above as a function of redshift for different chirp masses. The Poissonian 
constraint is plotted together with the expected rates in Figures |4] and [5] Since there are few to no close 
SMBH binary s ystems detected near z = 0, it can be assumed that e ^ 1. The only free parameter 
remaining in the Ijaffe & Backer! (|2003h model is the evolution index which determines the SMBH binary 
merger rate as a function of redshift. The grey regions in Figures [2]-|5] give the range of expected merger 
rates for — 1 < 7 < 3. The maximum and minimum rates found using the lSesana etaP J2008L l2009i) models 
are shown as thick dashed and thin dashed lines, respectively. Note, since the minimium predicted rate 
for Mc = 10^*^ Mq is not presented by the Sesana et al. models, this curve is not shown. Overall, the 
upper bounds obtained by pulsar timing data do not contrain the parameters of the SMBH binary merger 
models discussed in this paper beyond their currently accepted ranges. For the PPTA goal (20 pulsars timed 
at 100ns rms accuracy for five years), the results imply that either a detecti on will be rnade or 7 < 1-7 
at redshift 2: < 3. In order to pla ce constraints that will limit the models of ISesana et al.l (|2008l . l2009r) as 
well as the lJaffe & Backeii (|2003h model with 7 < — 1, one must either time 100 pulsars with 100 ns rms 



timing precision or 20 pulsars at the 10 ns level, both of which should be possible with the proposed Square 
Kilometre Array project. 

The Poissonian constraint is only useful for constraining the properties of the most massive SMBH 
binaries since these systems are rarer than their less massive counterparts and emit a stronger GW signal. 
Figure m shows that the ideal PPTA extended to 10 years of observations will just be able to place useful 
limits for the most massive systems if 7 were in the larger end of its possible range. It will be more interesting 
when we can time pulsars to the 10 ns level. Here, the Poissonian cons trai nt will be well be low that expected 
for ~10'^ Mq SMBH binaries from both the lSesana etaD(l2008[ 120091) and lwen et aD(l200^) models. Hence, 
with 20 pulsars timed with 10 ns rms timing accuracy for 10 years, we have a very good chance of detecting 
an individual source or will place very stringent constraints on r nodels of SMBH bi nary formation and 
evolution. These conclusions are consistent with the recent work of lSesana et al.l ( 2009 ). 
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Fig. 2. — Upper limits on the SMBH binary merger rate determined by the stochastic constraint discussed 
in §2.1 for different data sets: real data published by Uenet et al.l (120060 (open triangle), simulated data for 
20 PSRs-500ns-10yr (open square), 20 PSRs-100ns-5yr (cross), and 20 PSRs-lOOns-lOyr (open circle). 
The error bar is plotted as a dotted line when the constraint is in yalid. The filled gray area represents the 
expected region for the coalescence rate using the framework of IJaffe & Backer! (120031) together with the 
data from the SDSS (.Wen et al.,,2009i) with an evolution index — 1 < 7 < 3. The dashed l ines between 
< Igfl + z) < 0.7 indicate the maximum (thick) and minimum (thin) predicted rates from Isesana et all 



(2008, 



2009 ). 
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Fig. 3.— Same as Figured but for simulated data sets: 100 PSRs-lOOns-5 yr (filled triangle), 100 PSRs- 
lOOns-lOyr (filled square), 20 PSRs-lOns-lOyr (star), and 100 PSRs-lOns-10 yr (filled circle). 
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Fig. 4. — Upper limits on the S MBH merger rate using the Poissonian constraint discussed in §2.2 for 
different data sets: real data from benet et al. ( 2006 ) (open triangle), 20 PSRs-500ns-10yr (open square), 
20 PSRs-lOOns-5 yr (cross), and 20 PSRs-lOOns-lOyr (open cir cle). The filled gray a rea represents the 
expected region for t he coalescence r ate using the framework of IJaffe & Backen (120031) together with the 
data from the SDSS dWen et al.l 120091) with an evolution index, — 1 < 7 < 3. The dashed lines between 



< Igfl + z) < 0.7 indicate the maximum (thick) and minimum (thin) predicted rates from lSesana et al. 



(12008 



200% . Note that no upper limits on the coalescence rate are obtainable for SMBH of Mc = 10^ M( 







at these sensitivities. 
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Fig. 5. — Same as FigureH] but for the following simulated data sets: 100 PSRs-100 ns-5 yr (filled triangle), 
100 PSRs-100 ns-lOyr (filled square), 20 PSRs-lOns-lOyr (star), and 100 PSRs-lOns-lOyr (filled circle). 
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4. Summary 



We have shown that pulsar timing observations may be used to place constraints on the rate of coa- 
lescence of binary supermassive black holes distributed throughout the Universe. Two types of constraints 
were considered: a stochastic constraint and a Poissonian constraint. The stochastic constraint, which is 
based on a detection algorithm for the stochastic GW background, gives lower rates but it is only valid when 
the expected amplitude for a individual source is much less than the minimum detectable amplitude. When 
this is not the case, the Poissonian constraint must be used. This constraint is based on a continuous-wave 
detection algorithm and it assumes that the number of coalescence events are distributed according to a 
Poisson distribution. In both cases, it is assumed that the differential rate of coalescence varies sufficiently 
slowly over a range of chirp masses and redshifts. The precise numerical value of the constraint depends on 
the size of the interval over which the rate is assumed to be nearly constant. 

The implied rate constraint obtained from recently published data together with rate constraints ex- 
pected from future possible observing scenarios were compared to theoretical rates calculated from different 
models. It was shown that 20 pulsars timed with an accuracy of 100 ns, the go al of the PPTA project, will 



Wen et al. ( 


2009). The upt 


models dSesana et al. 


20od. 



(120091). The upper end of the range of backgrounds produced by SMBH population synthesis 



will be needed to further constrain these models if the GW background is not detected. It was also shown 
that if future observations can time a pulsar with 10 ns accuracy, a direct detection of one or more individual 
sources is highly likely. 
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